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The pressure and viscosity in two-dimensional sheared granular assemblies are investi- 
gated numerically for varying disks' toughness, degree of polydispersity and coefficient of 
normal restitution. 

In the rigid, elastic limit of monodisperse systems, the viscosity is approximately inverse 
proportional to the area fraction difference from (\> n ~ 0.7, but the pressure is still finite 
at 4> v . On the other hand, in moderately soft, dissipative and polydisperse systems, we 
confirm the recent theoretical prediction that both scaled pressure (divided by the kinetic 
temperature T) and scaled viscosity (divided by VT) diverge at the same density, i.e., the 
jamming transition point <f>j > cj} n , with the critical exponents —2 and —3, respectively. 
Furthermore, we observe that the critical region of the jamming transition disappears as the 
restitution coefficient approaches unity, i.e. for vanishing dissipation. 

In order to understand the conflict between these two different predictions on the diver- 
gence of the pressure and viscosity, the transition from soft to near-rigid particles is studied 
in detail and the dimensionless control parameters are defined as ratios of various time-scales. 
We introduce a dimensionless number, i.e. the ratio of dissipation rate and shear rate, that 
can identify the crossover from the scaling of very hard, i.e. rigid disks, in the collisional 
regime, to the scaling in the soft, jamming regime with multiple contacts. 



§1. Introduction 

One of the reasons for the growing interest in granular materials, i.e. collections 
of interacting macroscopic particles^ ^ ® ® ^ ® & ® -E) Mi Mi Mi Mi Mi Mi Mi Mi Mi Mi Mi Mi ,123 Mi Mi 
is the fact that these materials are different from ordinary matter.^ The pertinent 
differences do not preclude a description of (up to) moderately dense and nearly 
elastic granular flows by hydrodynamic equations with constitutive relations derived 
using kinetic theory EJ ,[HJ ,H7]( ,[28j) ,[29ji ,[30j ,[3TJ) ,[32j ,[33ji ,[3U) -\yhen n0 ntrivial correlations, 

such as long-time tails and long-range correlations, are present, one can apply fluc- 
tuating hydrodynamic descriptions to granular fluids and the latter can be obtained 
from kinetic theory as we \\m Mi Mi Mi Mi Mi Mi Mi Mi Mi 

Similar analysis cannot be applied to systems near the jamming transition. In- 
deed, we know many examples when the behavior of very dense flows cannot be 
understood by Boltzamnn-Enskog theoi^^^^^^^^ due to effects like 
ordering or crystallization, excluded volume, anisotropy and higher order correla- 
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tions. Therefore, to understand the rheology of dense granular flows, such as the 
frictional flowP and the jamming transition itself J^D an alternative approach is 
called for. 

Recently, Otsuki and Hayakawa have proposed a mean-field theory to describe 
the scaling behavior close to the jamming transition^ at density (area fraction) 
<fij. They predicted that both pressure and viscosity are proportional to (4>j — </>)~ 4 . 
Therefore, the scaled pressure, divided by the kinetic granular temperature T oc 
(4>j — (/>)~ 2 , is proportional to (cf>j — (j))~ 2 , while the scaled viscosity, divided by 
VT oc ((f>j — (p)^ 1 , is proportional to (4>j — 4>)~ 3 , irrespective of the spatial dimension. 
The validity of this prediction has been confirmed by extensive molecular dynamics 
simulations with soft disks. 

However, one can note that this prediction differs from other results on the di- 
vergence of the transport coefficients.^'^'^*''^''^' In particular, Garcia-Rojo et 
alP^ concluded that the viscosity for two-dimensional monodisperse rigid-disks is 
proportional to {(f) v — </>) -1 , where tp^ is the area fraction of the 2D order-disorder 
transition point, while the pressure diverges at a much higher <pp with p oc ((ftp — 
^-ljUlJiZJJiDJlSJlD ]y ot Qn iy i g i oca ti on Q f the divergence different, but also 
the power law differs from the mean field prediction in Refs. [39j) . l40|) . How can we 
understand these different predictions? One of the key points is that the situations 
considered are different from each other. As stated above, Garcia-Rojo et alP^'^J 
used two-dimensional monodisperse rigid-disks without or with very weak dissipa- 
tion, whereas Otsuki and HayakawsP^'^S discussed sheared polydisperse granular 
particles with a soft-core potential and rather strong dissipation. 

In order to obtain an unified description on the critical behavior of the viscosity 
and the pressure in granular rheology, we numerically investigate sheared and weakly 
inelastic soft disks for both the monodisperse and the polydisperse particle size- 
distributions. The organization of this paper is as follows: In the next section, 
we summarize the previous estimates for the pressure and the viscosity for dense 
two-dimensional disk systems. In Sec. El we present our numerical results for soft 
inelastic disks under shear in three subsections: In Sec. 13. 1\ the numerical model is 
introduced, Sec. 13.21 is devoted to results on monodisperse systems, and Sec. 13.31 to 
polydisperse systems. In Sec. 13.41 a criterion for the ranges of validity of the different 
predictions about the divergence of the viscosity and the pressure is discussed. We 
will summarize our results and conclude in Sec. HI 

§2. Pressure and viscosity overview 

In this section, we briefly summarize previous results on the behavior of pressure 
and viscosity in two-dimensional disks systems. Following Ref. 23J, we introduce the 
non-dimensional pressure 

P* = P/(nT) - 1, (2-1) 

where P is the pressure, n is the number density, and T = (m(v — (v)) 2 ) / (2N) is 
the kinetic temperature (twice the fluctuation kinetic energy per particle per degree 
of freedom) which is proportional to the square of the velocity fluctuations of each 
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particle. We also introduce the non-dimensional viscosity 

rf = n/{pv T s /2) (2-2) 

where p denotes the particles' material density, p B = p<fi is the bulk area density, 
the fluctuation velocity is denoted by vt = W2T/m, so = \/27rcr/8, the mass of a 
grain (we assume all grains to have the same mass) is denoted by m, and the mean 
diameter of a grain (disk) is denoted by a. It should be noted that pvxSo/2 is the 
viscosity for a monodisperse rigid-disk system in the low-density limit and correct to 
leading order in the Sonine polynomial expansion. For later use, we also introduce the 
mean free time tE which is defined as the time interval between successive collisions. 
This leads to the collision rate t E l = VT<t>g{4>) / 'sq = vt/X in the case of dilute and 
moderately dense systems of rigid disks, where A is proportional to the mean free 
path. 

In the first part of this section, let us summarize previous results for elastically 
interacting rigid disk systems. In the second part of this section, we show other 
previous results for soft granular disk systems under shear. 

2.1. Rigid disk system in the elastic limit 

For the equilibrium monodisperse rigid-disk systems, the reduced pressure P* 
of elastic systems at moderate densities cj) < 0.67 is well described by the classical 
Enskog theoryS»H3Jira 

PI = 2&fe(0). (2-3) 
with the aid of improved pair-correlation function at contact 

</> 3 /16 

=&»(*)- 8( 1 ■ ( 2 ' 4 ) 

where g2(<ft) = T[~^pr m Eq. (|2-4p was proposed by Henderson in 1975.^ In the 
regime of high density <f> > 0.65, the reduced pressure becomes, first, lower than 
(|2-3j) because of ordering (crystallization) and, second, diverges at a density 4>p due 
to excluded volume effects. This behavior is quantitatively fitted by 

^dense = ^~iK^P ~ <t>) ~ ^ (2'5) 
<PP - <P 

with 4>p = it/ (2y3), h(x) = 1 + c\x + C3X 3 , and the fitting parameters c\ = —0.04, 
and c 3 = 3.25.EU'Si'S3'153 As shown in references ,HBJ) , ED an interpolation law 
between the predictions for the low and the high density regions: 

Pq = PI + MmPl nse - PH (2-6) 

with M(cf)) = [1 + exp(-(0 - <£ c )/m ] -1 , <p c = 0.699, and m = 0.0111, fits well the 
numerical data for P* . The quality of the empirical pressure function Pq is perfect, 
except for the transition region, for which deviations of order of 1% are observed in 
the monodisperse, elastically interacting rigid disk system. 
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The dimensionless viscosity for monodisperse elastically colliding rigid disks is 
well described by the Enskog-Boltzmann equation 



Ve 



' + 20+ (l + -U 2 <7 2 (</>) 



92(9) V ^ 



(2-7) 



Note that g2(4>) satisfies 52 (<^) ~ 24(0) ~ 9q{4>) = Pq I {^9)1 f° r 9 ^ <?V A dominant 
correction, see Eq. (|2-8|) below, controls the viscosity for higher densities, closer to 



9 ~ 9v 

Equation (|2-7|) can be used for low and moderate densities, but it is not appro- 
priate close to the crystallization area fraction ^ c j3DJ3UJ33J31JilS)J3n)JlSJ Therefore, 
an empirical formula for 77* has been proposed as 

= f 1 + A ~ c f) 1* (2 ' 8) 

V 9^~9 9vJ 

which can fit the numerical data for < <j) < 4> v with two fitting parameters c v = 
0.037 and 9^ = 0.7lP^ Note that the last term is an improvement of the original 
empirical frP» that makes rf L approach unity for — > 0. Note that 77* in Ref. [36]) 
was obtained from a non-sheared system by using Einstein-Helfand relation.^ 

A slightly different empirical form for the non-dimensional viscosity was pro- 
posed by KhairPSJ (based on simulations of a sheared system): 



V*K= l + T^r f )V%, (2-9) 




with the same and (j)^ as before. The reasons for the difference between the 
viscosity in a sheared and a non-sheared system is an open issue and will not be 
discussed here. 

We also introduce the scaled temperature given by 

Til - e 2 ) 

T* = { „/ (2-10) 

m / j z SQ 

for sheared inelastic rigid-disks, where e and 7 are the coefficient of restitution and 
shear rate, respectively. Luding observed that the empirical expression 

T* K = T^TT (2-11) 

fits best the numerical data for monodisperse rigid disksp^ while 

T* L = (2-12) 
9 92{9) 

slightly overpredicts the scaled temperature. 

For polydisperse elastic rigid-disk systems, many empirical expressions for the 
reduced pressure P* have been proposed, see e .gpBES)EDEIO>'El> ^ s known that P* 
diverges around max — 0.85 for bi- and polydisperse systems, but there is no theory 
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to our knowledge that predicts the dependence of <^ max on the width of the size 
distribution function that was observed in rigid-disk simulations PSJ >S9 Dependent 
on the dynamics (rate of compression) , on the material parameters (dissipation and 
friction), and on the size-distribution, different values of <^ max can be observed. In 
several studies, the critical behavior was well described asymptotically by a power 
law 

P*d ~ (<Amax " (P)' 1 (2-13) 

see Refe. 1231,113), HSU . 

No good empirical equation for the viscosity of polydisperse rigid-disk systems 
in the elastic limit has been proposed to our knowledge. However, if we assume 
that the viscosity behaves like that of the monodisperse rigid-disk system, we can 
introduce the empirical expression 

V*d ~ (0max " 0)^ (2-14) 

as a guess. Here, we assume that the pressure P* and the viscosity rj* for the 
polydisperse system diverge at the same point (/> max , which differs from the case of 
the monodisperse system, where P* and rf diverge at different points 4>p and <f>^ 
due to the ordering effect. 

2.2. Soft- disk system 

Let us consider a sheared system of inelastic soft-disks characterized by the 
non- linear normal repulsive contact force k8 with power A, where k and 5 are the 
stiffness constant and the compression length (overlap), respectively. For this case, 
Otsuki and Hayakawa-^'HS proposed scaling relations for the kinetic temperature 
T, shear stress S, and pressure P, near the jamming transition point cftj ~ 0.85: 

T = l^ r ± (i^) ■ s = W v * s ± (^) > p = W y '* v ± (jijs) > ( 2a5 ) 

where <P = <fi — is the density difference from the jamming point. This scaling 
ansatz is based on the idea that the system has only one relevant time-scale r ~ 
diverging near the transition point cj)j, and the behavior of the system is dominated 
by the ratio of the time scale r and the inverse of the shear rate j. This idea is often 
used in the analysis of critical phenomena. 

The scaling functions T+(x), S + (x), and V+(x) satisfy 

lim 7~+(x) = x, lim S+(x) = 1, lim V + (x) = 1 (2'16) 

for (J) > <j)j, i.e., for higher area fraction. The pressure and shear stress scaling - in 
this limit - represent the existence of a (constant) yield stress S = Sy- The scaling 
for the temperature is obtained from the assumption that a characteristic frequency, 
oj = jS/i^nT), is finite when 7 — > in the jammed state <j> > <pj, see Ref. [57]) . E3 

*' Here, we should note that oj is proportional to the Enskog collision rate ui = (1 — e 2 )t~j s 1 /2, see 
Ref. I24p . in the unjammed state well below the jamming point, <f> < 4>j, i.e., in the collisional flow 
regime. Due to the prefactor (f — e 2 )/2, we can identify u> with the characteristic dissipation rate. 
The different time-scales (inverse frequencies) and their relative importance are discussed below in 
subsection 13. 1.21 
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On the other hand, for lower area fraction, T-(x), S-(x), and V—(x) satisfy 

lim T-{x) = x 2 , lim S-(x) = x 2 , lim V-(x) = x 2 (2-17) 

for (j) <C 4>j, which represent Bagnold's scaling law in the liquid phase. 

Furthermore, for diverging argument x, i.e., at the jamming point J with <P — > 0, 
the scaling functions T±(x), S±(x), and V±(x) should be independent of ^ and thus 
satisfy: 

lim T±(x) = x x * /a , lim S±(x) = x y * /a , lim V±(x) = x y '* /a . (2-18) 

x— >oo x— >oo x— >oo 

The critical exponents in Eq. ([2- 15j) are given by 

A + 4 

x$ = 2 + A, y$ = y' <p = A, and a = — - — , ( 2 '19) 

which depend on some additional assumptions,® such as the requirement that the 
pressure P for 7 — > 0, in the jammed state ^ > 0, scales with the force power-law as 
P~# 4 , see Refe-ESIlJEfl!. 

Thus, the temperature T, the shear stress S, and the pressure P, below the 
jamming transition point in the zero shear limit 7—7-0 obey: 

T ~ {4>j — 0)~ 2 7 2 , S ~ {<t>j - 0)" 4 7 2 , P~(^j-0)- 4 7 2 - (2-20) 

Both the viscosity 77 = 5/7 and pressure P, at the jamming transition point, diverge 
proportional to the area fraction difference to the power —4. Substituting Eqs. (|2-20p 
into Eqs. (|2-ip and (|2-2p . the reduced pressure P* and the dimensionless viscosity 
rj*, in the vicinity of the jamming point are respectively given by 

pj ~ - 4>r 2 , (2-2i) 

r/} ~ (0J - <P)- 3 . (2-22) 

It is remarkable that the scaling relations (|2-20p ^ (|2-22p below the jamming tran- 
sition point are independent of A, even though the exponents in Eq. (|2-19p depend 
on A. The validity of Eqs. (|2-2ip and (|2-22|) for various A has been numerically ver- 
ified.®'^ However, the conjecture that the scaling relations (|2-2ip and (|2-22p are 
applicable in the hard disk limit seems to be in conflict with the empirical relations 
Eqs. (|2-13|) and (|2-14|) for elastic rigid-disk systems. 



§3. Numerical results 



In this section, we numerically investigate the reduced pressure P* and viscosity 
77* of sheared systems with soft granular particles, with special focus on the rigid- 
disk limit. In the first part, our soft-disk model is introduced. In the second part, 
we present numerical results for monodisperse systems, while in the third part the 
results for polydisperse systems are presented. 
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3.1. The soft-disk model system 

3.1.1. Contact forces and boundary conditions 

Let us consider two-dimensional granular assemblies under a uniform shear with 
shear rate 7. Throughout this paper, we assume that granular particles are fric- 
tionless, without any tangential contact force acting between grains. For the sake of 
simplicity, we restrict ourselves to the linear contact model with A = 1. We assume 
that all particles have identical mass regardless of their diameters. The linear elastic 
repulsive normal force between the grains i and j, located at r, and rj, is: 

fe\{rij) = kO (aij - rij) (o-ij - rij), (3-1) 

where k and r^ are the elastic constant and the distance between the grains r^ = 
\ r ij\ = \ r i — r j\i respectively. oij = (<7j + o~j)/2 is the average of the diameters of 
grains i and j. The Heaviside step function 0{x) satisfies 0(x) = 1 for x > and 
0(x) = otherwise. The viscous contact normal force is assumed as 

where £ is the viscous parameter. Here, «y )n is the relative normal velocity between 
the contacting grains Vij,n — (vi — ^j) * l^ij 5 where Vi and Vj are the velocities of 
the centers of the grains i and j, respectively. 

In order to obtain a uniform velocity gradient 7 in y direction and macroscopic 
velocity only in x direction, we adopt the Lees-Edwards boundary conditions. The 
average velocity c(r) at position r is given by c(t) — ^y€, X) where €x,a. 

is a unit 

vector component given by e x>a = S xa , where a is the Cartesiani coordinate. 

3.1.2. Discussion of dimensionless quantities 

There are several non-dimensional parameters in our system. One is the resti- 
tution coefficient e given by 



cxp 



y/2k/m - (C/m) 2 <>XPI U ' ] ' ' :> " ! 

with the pair-collision contact duration t c = -k/ \j2kjm — (Q/m) 2 . Another is the 
dimensionless contact duration 

r* = t c7 (3-4) 



that represents the ratio of the two "external" time-scales of the system ^2 • "Exter 



nal" means here that these time scales are externally controllable, i.e., the contact- 



*^ The contact duration t c is well denned for two masses connected by a linear spring-dashpot 
system and corresponds to their half-period of oscillation. A particle in a dense packing (connected 
to several masses by linear spring-dashpots) has a somewhat higher oscillation frequency, but the 
order of magnitude remains the same. Particles with non-linear contact models can have a pressure 
dependent t c , but are not considered here. 

**' One can see r* = (a A /)/(a/t c ) also as the ratio of the two relevant velocities in the dense 
limit, i.e., as the ratio of the local velocity of horizontal layers that are a diameter of a grain, a, 
apart, and the local information propagation speed a/t c in a dense packing. However, the ratio 
of velocities makes only sense in the dense, soft regime, since t c is not a relevant time-scale in the 
dilute, near-rigid regime. 
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duration is a material parameter and the inverse shear rate is externally adjustable. 

In all cases studied later, we have r* <C 1, which means that the shear time scale 
is typically much larger than the contact duration, i.e., we do not consider the case 
of very soft particles, which is equivalent to extremely high shear rates. Therefore, 
r* will be used as dimensionless control parameter in order to specify the magnitude 
of stiffness: The rigid disks are reached in the limit r* — > 0. 

The third time-scale, £g, in the system is an "internal" variable, i.e., cannot 
be controlled directly. This time scale is proportional to the inverse characteristic 
frequency of interactions, i.e., the mean free time, tg, in the dilute case or the rigid- 
disk limit. This defines the (second) dimensionless ratio of times 

T* E = t E j (3-5) 

relevant in the dilute, collisional regime. 

The third dimensionless number is defined as the ratio of contact duration and 
mean free time, 

= £ = (3-6) 

see Eq. (53) in Ref. |2^|) . The meaning of this dimensionless number is as follows: 
For very small t* e <C 1 one is in the binary collision regime, for large t* e > 1, one 
is in the solid-like regime with long-lasting multi-particle contacts. In the hard disk 
limit r* —?• 0, we can identify t* e with the coordination number as will be shown in 
Fig. [8l Namely, finite r* E in the near-rigid situation means that the system is in a 
jammed phase. 

The binary collision regime, t* e — > 0, cannot be controlled directly, since t E 
is a function of temperature, which depends on e and 7. On the other hand, the 
rigid-disk limit, r* — > 0, can be approached/realized by either (i) vanishing shear 
rate, 7 — > 0, or (ii) near-rigid particles with high stiffness, k — > 00 (with controlling 
the variable £ to maintain a constant restitution coefficient e). 





ratio of times 


ratio of velocities / stresses 


regime of relevance 


r* 




Vaj/vc = aj/(a/t c ) 


near-rigid, high density (a 3> A, t c 3> tE) 


T E 


w-r 1 

tolls 


v\j/vE = X-y/(X/t E ) 
ve/vc = a/t E /(cr/t c ) 


rigid, low density (a <C A, t c <g ts) 
near-rigid, low and moderate densities 


t uj 


tc/uj- 1 


nT/S = 2T E /(l~e 2 ) 
t c jS/(nT) = t;/t* 


well denned in sheared systems 
well defined in all systems 



Table I. Summary of the dimensionless numbers discussed in the text, where i c , 7 , tE, are 
contact duration, inverse shear rate, mean free time, and inverse characteristic dissipation rate, 
respectively. The velocities VL-y, v c , and ve are the shear velocity of layers separated by length 
L, the speed of sound propagation in a dense packing, and the speed of sound propagation in a 
dilute packing, respectively. The relevant lengths L can be the diameter a (in the dense limit), 
the mean free path A = A(</>) (in the dilute limit), or their sum (for all densities). 



Furthermore, we can introduce dimensionless numbers that are related to the 
inverse characteristic dissipation rate w _1 EE which has the meaning of the energy 

*' Note that the identity a; -1 = 2^/(1 — e 2 ) is true in the dilute, collisional limit only. For 
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dissipation time-scale. For e — > 1, dissipation is becoming very slow, while for small 
e ~ 0, considerable energy can be dissipated, within a time of order of ts or t c . 
Replacing by u^ 1 in Eqs. (|3-5|) and (|3-6|) . we obtain 

r^7/«, (3-7) 
= ^ . (3-8) 

It should be noted that r* and r*^ approximately satisfy the relations r* ~ 2te/(1 — 
e 2 ) and r* w « (1— e 2 )r* E /2, respectively, in the collisional regime, where the prefactor 
plays an important role, as will be demonstrated later. 

The consequences of the interplay among these dimensionless numbers will be 
clarified and discussed in the following sections. Furthermore, we will identify the 
dimensionless number that - we believe - allows us to distinguish between the two 
scaling regimes. 

3.1.3. Simulation parameters 

We examine two systems with different grain diameters and composition. The 
first monodisperse system consists of only one type of particles, whose diameters are 
o~q. The other polydisperse system consists of two types of grains, and the diameters 
of grains are 0.5<7o, and do, where the numbers of each type of grains are 0.8N and 
0.2N, respectively, with the total number of particles N. The reasons to study such 
a polydisperse system are (i) to avoid crystallization and (ii) to compare our new 
near-rigid data with previous results from rigid disks.^'^ 

In our simulations, the number of particles is N = 2401 except for the data in 
Figs. [11] and [T2l where we have used N = 20000. We use the leap-frog algorithm, 
which is second-order accurate in time, with the time interval At = 0.2ym/fc. We 
checked that the simulation converges well by comparison with a shorter time-step 
At = 0.02^/mJk. 

The pressure and the viscosity are respectively given by 

N N 



p = \ E E r v L/eiM + Usinj, %>)] + E / ' (3 ' 9) 

N \ 

[Unj) + f vis (rij, n)] + £ ) , (3-10) 




8=1 / 

with the volume of the system V, the relative distance vector = (rij tX , rij iV ), with 
T%j = and the peculiar momentum •p i = (pi iX , pi jV ) = m(vi — A fyie x ). 

3.2. Mono-disperse system 

In Figs. ID^a) and (b), we plot P* as a function of the area fraction <j) in the 
monodisperse system with e = 0.999 for < (j) < 0.6 and 0.5 < <fi < 0.9, respectively. 
Most of all data of P* seem to converge in the rigid-disk limit (r* — > 0). Moreover, 



higher densities and for softer particles, one hasw" 1 > 2t E /(l-e 2 ), i.e., energy dissipation becomes 
somewhat slower when approaching the jamming transition. This is consistent with a slower energy 
decay due to the reduced dissipation rate, proposed in Eq. (52) in Ref. 12 4[) 
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the data for P* with < 0.6 are consistent with Pq, see Fig. H^a), while P* for 
4> > 0.7 in Fig. QJb) deviates from Pq in the soft case of r* = 1.11 x 10~ 3 , and 
also in the rigid-disk limit. Only the simulations with r* = 1.11 x 10 -4 are close to 
Pq - seemingly by accident. At high densities, for very soft particles, the stress is 
considerably smaller than predicted by Pq, while for near-rigid particles, we observe 
a higher stress. 



(a) 





Fig. 1. The reduced pressure P* as a function of the area fraction <j> in the monodisperse system 
with e = 0.999, for different r*, as given in the inset, and for cfi < 0.6 (a) and (f> > 0.5 (b). 




In order to check the possibility that the restitution coefficient is the reason for 
the deviation between the numerical data and PX in Fig.[T^b), we plot P* for different 
e, for t* = 1.11 x 10~ 5 in Fig. [2j At high densities, for inelastically interacting 
particles, e = 0.99, the stress is considerably smaller than predicted by Pq, while 
for more elastic particles, we observe a higher stress. Only the almost elastic case 
e = 0.9999 is close to the prediction. 

The low pressure for e = 0.99 is due to the existence of a shear-band - see below. 
For all other situations, no shear-band is observed, however, different patterns of 
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defect lines in the crystal are evidenced for e = 0.9990 and e = 0.9995, while an 
almost perfect crystal is observed for e = 0.9999, where slip-lines appear. It should 
be noted that the positions of the slip-lines (shear-bands of width W = d) don't 
move in the steady state of one sample, but vary among different samples. 

(a) (b) 



0.02 ■ 

U' 

0.01 " 

■ 

^^^^^^^^^^^^^^^^ 

■ 

■ 

0.01 ■ 
-0.02 

-20 20 

y' 





Fig. 4. (a) The scaled velocity u' (like in Fig. 0J for (p = 0.84, r* = 1.11 x 10" 5 , and e = 0.999. 
(b) Snapshot of the monodisperse system from (a). 

We have confirmed the existence of shear-bands for (p > 0.7 with e = 0.99 in 
Fig. [31 We plot the velocity u(y) in x direction as a function of y for <fi = 0.76, 
r* = 1.11 x 10 -5 , and e = 0.99 in Fig. [3] (a), where the velocity gradient exists only 
in the regions y/o~o < —20 or y/o~Q > 20. The apparent inhomogeneity is observed 
in the snapshot of the system, see Fig. E^b). On the other hand, such a shear-band 
could not be observed for the case of e = 0.999. Note that the shear-band formation 
in our system is different from that for the dilute case®^ in which dense strips align 
at 45 degrees relative to the streamwise direction. Fig. [3] shows that the system is in 
an uniformly sheared state with some density fluctuations, see Fig. Iljb). Actually, 
here deformations take place irregularly and localized - together with defects and 
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Fig. 5. (a) The scaled velocity it' (like in Fig. © for <p = 0.84, r c * = 1.11 x 1CT 5 , and e = 0.9999. 
(b) Snapshot of the monodisperse system from (a). 



slip planes - so that the velocity profile looks smooth and linear only after long-time 
(or ensemble) averaging. For the case of e = 0.9999, almost perfect crystallization is 
observed, but slip-lines exist, see Figs. G^a) and (b). 

This is in conflict with the observations of Ref. 24), where shear-bands were 
observed at densities around (ft « 0.70, (ft « 0.73, and (ft « 0.78, for e > 0.99, 
e = 0.95, and e = 0.90, respectively. In this paper, for the case of e = 0.999, 
no shear band is observed, however, in the simulation of the sheared inelastically 
interacting rigid-disks with e = 0.998 in Ref. 24), a shear band was reported. 

We identify two differences between the systems in this paper and Ref. I24p . The 
first difference is the softness of the disks that, however, should not affect the results 
as long as we are close to the rigid-disk limit. The second difference is the protocol 
to obtain a sheared steady state with density (ft. In this paper, first an equilibrium 
state with density (ft is prepared and then shear flow and dissipation between the 
particles is switched on to obtain the sheared steady state. In contrast, in Ref. 
|24"|) . the system of sheared inelastically interacting disks was studied by slowly but 
continuously increasing the density (ft. 

The dimensionless viscosity rf for monodisperse systems with e = 0.999, and 
different r* is shown in Fig. We note that both P* and r/* converge for more 
rigid disks r* — > 0, but not to the empirical expression rf L from Eq. (|2-8p . It can be 
used in a wide range of (ft, as one can see in Fig. [5(b), but - even though behaving 
qualitatively similar - the numerical data clearly deviate from rf L : For (ft > 0.7, 
in the rigid-disk case, rf L diverges at (ft^ = 0.71, whereas rf in the near-rigid case 
exponentially grows like the Vogel-Fulcher law, which remains finite above (ft^. 

The difference between the numerical data for rf and rf L results from both 
elasticity and dissipation, as shown in Fig. [71 where the dependence of rf on (ft for 
r* = 1.11 x 10~ 5 and different coefficients of restitution e are plotted. The viscosity 
rj* , like the pressure P*, approach rf L and Pq in the elastic limit e — > 1, i.e., they 
converge to the results of the elastic rigid-disk system. It should be noted that 
Figs. E)^a) and [7(a) suggest that the singularity around (ft = (ft n is an upper limit, 
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Fig. 6. (a) The dimensionless viscosity rf as a function of the area fraction <j> in the monodisperse 
system for e = 0.999, and different r* as given in the inset, (b) rf jr\* E as a function of the area 
fraction from the same simulations as in (a). 



only realized in the rigid disk limit and for e — > 1. As will be discussed below, 
for given r* and e, the simulations deviate more and more from the rigid disk case 
with increasing density. The smaller r*, i.e., the stiffer the disks, the better is the 
upper limit approached - but for finite dissipation and for near-rigid disks, there is 
always a finite density where the elasticity (softness) becomes relevant and leads to 
deviations from the upper limit. Above that density, it seems that the divergence 
of the viscosity takes place at the same point as the pressure, and another inverse 
power law can be a fitting function for (j) < <p„. 
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Fig. 7. The dimensionless viscosity rj* as a function of the area fraction (j) in the monodisperse 
system for r c * = 1.11 x 10" 5 and e = 0.99, 0.999, 0.9999. (b) rj* /rj* E as a function of the area 
fraction from (a). 



In rigid-disk systems, the coordination number Z should be identical to zero 
because the contacts between the particles are instantaneous. Hence, in the rigid- 
disk limit of soft-disks, it is expected that the coordination number Z vanishes, which 
is confirmed by Fig. E{a). Here, it should be noted that the coordination number 
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Z is almost identical to the dimensionless number t* e W^ Indeed the relationship 
Z « t* e can be verified in Fig. |H](b) , where we plot the ratio t* e /Z as function of 
the area fraction <j> for monodisperse systems with e = 0.999 and several r*. Here, 
we have measured the coordination number as 



Z = ^ Z ^{e(a ij -r l3 ))/N . 

i j^i 

If we use the mean-field picture, we can understand the relation Z 
in Appendix [XI 



(3-H) 



t* e as shown 
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Fig. 8. (a) The coordination number Z plotted as function of the area fraction <f> for monodisperse 
systems with e = 0.999 for several t* values, (b) t* e /Z plotted as function of the area fraction 
4> from the same simulations as in (a). 

We also show the scaled temperature T* for the soft-sphere monodisperse system 
in Fig. As expected, T* approaches the empirical expression in Eq. (|2-11|) . 
This result also supports our conjecture that the rigid-disk limit of the soft-disk 
assemblies coincides with the rigid-disk system when the coefficient of restitution e 
is sufficiently close to unity. 

3.3. Poly-disperse systems 

In order to understand the polydisperse situation, we also study systems with 
different r* and different e values - as in the previous subsection. The reduced 
pressure P* and the dimensionless viscosity rf are almost independent of r* and e 
for moderate densities (</> < 0.8), as shown in Fig. [TUl where P* and 77* are plotted 
as functions of the area fraction <f>. For low densities, the simulation results of P* 
agree with the scaling given by PI, while the asymptotic scaling behavior of rf is 
described by rj^ only above eft ~ 0.8. Here, we have used (p max = 0.841 for P* and rj* 
in Eqs. flZZED and 

However, when looking more closely, there are distinct differences between P* 
and P£, and between 77* and n* d for <p > 0.83. In Fig. [T[J P* and rf are plotted 
from polydisperse systems with rather strong dissipation, e = 0.9, where we have 
used particular values for (/> max = 0.841 and <j)j = 0.8525 in order to visualize their 
different behavior. Although P* is still finite for <p > ^max m the hard disk limit, 
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Fig. 9. (a) The scaled temperature T* as a function of the area fraction <f> for the monodisperse 
system at = 1.11 x 10 -5 and different e. (b) T* /T% as a function of the area fraction </> from 
the same simulations as in (a). 
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Fig. 10. (a) The dimensionless pressure P* as a function of the area fraction (j> for polydisperse 
systems with several different r* and e, where we have used <^ max = 0.841 for PJ and rf d . The 
prefactor for PJ oc (</> max — <A) _1 is chosen as 2^> max P^'^'^ (b) The dimensionless viscosity 
rf as a function of the area fraction cj> from the same simulations as in (a). Here, we have used 
the prefactor 7.0 for 77^ oc (0 max — </>) _1 . 



even for the smallest t c values, both P^ and diverge at c/> ma x &s (^max 

6)- 1 . On 

the other hand, in the same high density range, P* and 77* are consistent with PJ 
([2^2TD and n* ([2^2l P'BB i n the rigid-disk limit (r* = 1.11 x 1(T 6 ), as will be shown 
below. 

In order to verify whether the critical behavior of P* and 77* can be described 
by Pj and 7/}, we plot P* and rj* as functions of cpj — (ft in Fig. [12l Here, we plot 
only the data for <p < <pj because we discuss the scaling behavior of P* and rf in 
the unjammed regime in this paper. P* and 77* in the rigid-disk limit approach PJ 
and rjj, which satisfy ((f) j — 4>)~ 2 and (4>j — 4>)~ 3 , respectively. 

It should be noted that the plateaus in Fig. [121 close to the jamming transition 
point, for eft ~ cpj, can also be predicted from the scaling theory, by rewriting Eqs. 
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Fig. 11. (a) The dimensionless pressure P* as a function of the area fraction 4> for the polydisperse 
system for e = 0.9 and several r*. (b) The dimensionless viscosity rf as a function of the 
area fraction <f> in the polydisperse system from the same simulations as those in (a). Here, we 
have used <f)j — 0.8525 for Pj and rfj, and <j) max = 0.841 for and 77^. The prefactors for 
P d * oc (0 m ax-0) _1 and r)* d oc (0 max -0) _1 are 2<^) max and 7.0, respectively. For P} oc (<j) j ~ <f))~ 2 
and 77} oc (<jij — </>)~ 3 > the prefactors are chosen as 0.07 and 0.002, respectively. 



(|2-15|> — 1|2- 19|> . More specifically, the arguments are taken to the power — 1/a: 

T = r* /a n (^L) , s = ^' a s' ± (w\ , p = ^' a v ± (w\ > (342) 



where we have introduced 7±{x) = x~ x,p T±(x~ a ), S'±(x) = x~ y,p S±(x~ a ), and 
V±(x) = x~ y '$V±(x~ a ). The scaling functions satisfy lim^ _>>o T±(x) = lim, r _ s .o S±(p) = 
lim x ^oV'±(x) = const. Substituting these relations into Eqs. (|2-ip (|2-2p . with Eqs. 
(|2-19p . 7/ = 5/7, Z\ = 1, and the definition of r* given by Eq. (|3-4p . the scaling 
relations of P* and n* are obtained as 

p * = rt-^Vl , = r^U\ . (3-13) 

Here, the scaling functions satisfy linx^o ^±0*0 = h m x->o %±( x ) = const. Therefore, 
the plateau for P* and r/* in Fig. □Jshould be proportional to (1/r*) 4 / 5 and (1/t*) 6/5 , 

*4/5 *6/5 

respectively, which is confirmed by Fig. [131 where we plot P*t c and tj*t c as a 

function of (4>j — <j>)/ t* 2 ^ ■ 

Whether the simulation pressure is described by PJ or PJ, and whether the 
viscosity is given by 77^ or 77}, strongly depends on the coefficient of restitution e. In 
Figs. [T4"ttl6| we plot P* and rf as functions of </> for various e, involving the very 
high dissipation case e = 0.1, an intermediate case e = 0.99, and a low dissipation 
case e = 0.998. Using fitting values </> max = 0.841, 0.848, and 0.851, based on a fit 
starting from very low densities, corresponding to various e = 0.1, 0.99 and 0.998, 
respectively, we can approximate the data of P* best by PJ = 20 max /((/> max — <p). 
On the other hand, we assume that <pj is independent of e, and fix <pj = 0.8525 for 
all e, as confirmed this by our numerical simulations. 
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Fig. 12. (a) The reduced pressure P* plotted as a lunction ol <j>j — <fi for polydisperse systems with 
e = 0.9 and several r* based on the simulations used for Fig. 1111 (b) The dimensionless viscosity 
rf from the same simulations as those in (a). Here, we have used <j)J = 0.8525 for P} and rfj, 
and ^ ma x = 0.841 for P^ and n%. 
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Fig. 13. (a) Plots of P*t* 4 ^ 5 versus (cj>j — <^)/r* 2//j for polydisperse systems with e = 0.9 and 
several r*. (b) Plots of r/*r* 6//5 versus (4>j — {/>)/r* 2//j for polydisperse systems with e = 0.9 and 
several t*. 



Even in the case of strong inelasticity (e = 0.1), as shown in Fig. [TU PJ and rjj 
characterize the behavior of P* and rf near the jamming transition point, while 
and rjj deviate for (j) > 0.83. The range where Pj and r/j characterize the pressure 
and the viscosity becomes narrower as e — > 1, while the range of validity of P^ 
becomes wider, as shown in Figs. [15] and [161 For e = 0.998 (Fig. [16]) . the difference 
between P^ and Pj appears only in a small region of 4> which is shown in Fig. [T71 

Since the scaling behaviors of P* and rf agree with Pj and n*j near (fij, we 
conclude that the critical behavior for inelastic near-rigid systems is well described 
by Pj and rfj , as proposed in Refs. [3"9"J) .FfO" |) . The scaling plot in Fig. Q2] supports 
the validity of the critical behaviors concerning both the plateaus and the lower 
densities. However, such predictions cannot be used for almost elastic and perfectly 
elastic systems, neither mono- or polydisperse, whose critical behavior is described 
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by P£ and n* d instead. 



(a) (b) 




0.80 0.83 0.86 0.80 0.83 0.86 



Fig. 14. (a) The reduced pressure P* as a function of the area fraction <j> for the polydisperse 
system with e = 0.1 and several r*. (b) The dimensionless viscosity rf from the same data as 
those in (a). We have used b = 0.07 and 4>j = 0.8525 for Pj and r\j, and max = 0.841 for Pd 
and r\ d . We used 4>.j = 0.8525 for P} and rfj, and max = 0.841 for P d and rj* d The prefactors 
for Pj oc (0 m ax — (^) ~ 1 and rf d oc (<£ max — are 20 max and 7.0, respectively. The prefactors 
for Pj oc (4>j — 0) -2 and rf oc (0j — 0) -3 are given by 0.07 and 0.002, respectively. 
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Fig. 15. (a) The reduced pressure P* as a function of the area fraction 4> f° r the polydisperse 
system with e = 0.99 and several r*. (b) The dimensionless viscosity rf obtained from the same 
data as those in (a). We have used <f)j = 0.8525 for Pj and rjj, and max = 0.848 for Pd and rjd- 
The prefactors for P2 oc (<^ max — an d Vd (0max — 'W -1 are 20 max and 10.0, respectively. 
The prefactors for P} oc (cf>j — <j))~ 2 and ??} oc ((f).j — 4>)~ 3 , are 0.035 and 0.0015, respectively. 



3.4. Dimensionless numbers and a criterion for the two scaling regimes 

In Sec. 13.31 we reported a crossover from the region satisfying Eqs. (|2-13p and 
(|2-14j) to the region satisfying Eqs. (|2-2ip and (|2-22|) . Figure [T51 presents a schematic 
phase diagram in the plane of the restitution coefficient e and the area fraction <p, 
where -1 denotes the region satisfying the scaling relations given by Eqs. (|2-13p and 
(|2-14|) . and OH denotes the region satisfying the scalings given by Eqs. (|2-21[) and 
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Fig. 16. (a) The reduced pressure P* as a function of the area fraction <f> for the poly disperse 
system for e = 0.998 and several r*. (b) The dimensionless viscosity 77* obtained from the same 
data as those in (a). We have used (f>.j = 0.8525 for P} and rjj, and m ax = 0.851 for P2 and rjd- 
The prefactors for P£ oc (<^ max — </ , ) _1 and 77^ oc (<£ max — <£) _1 are 20 max and 25.0, respectively. 
The prefactors for P} oc ((f>j — <j))~' 2 and rfj oc (</>./ — <?!>) _3 , are 0.01 and 0.001, respectively. 



(a) (b) 




Fig. 17. (a) The reduced pressure P* plotted as a function of cj)j — <j> for polydisperse systems 
with e = 0.998 and several r* based on the simulations used for Fig. 1161 (b) The dimensionless 
viscosity 77* obtained from the same simulations as those in (a). We have used cj>j — 0.8525 
for Pj and rjj, and (4 m « = 0.851 for P2 and rja- The prefactors for P^ oc (<^ max — and 
rf d oc (<^> m ax — are 20 max and 25.0, respectively. The prefactors for Pj oc (<fij — <^) 2 and 
f)*j oc (</>./ — <?i)~ 3 , are 0.01 and 0.001, respectively. 



(|2-22|) . For each e, the high density region satisfies Eqs. (|2-21|) and (|2-22|) . while 
the low density region satisfies the scalings given by Eqs. (|2-13p and (|2-14p . As the 
restitution coefficient approaches unity, the region of OH becomes "narrower", and 
disappears in the elastic limit. 

Now, let us discuss which of the dimensionless numbers t^, t* e , t* or r* w can 
be used as the criterion to distinguish between the two scaling regimes. It should be 
noted that the dimensionless number for the criterion must be a monotonic function 
of (f), because the scaling relations Eqs. (|2-13p and (|2-14p appear in the higher density 
region and the scaling relations Eqs. (|2-21|) and (|2-22[) appear in the lower density 
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jammed 




<|>j 

Fig. 18. A schematic phase diagram of the region (-1) satisfying Eqs. (|2-13p . (|2-14p and the region 
(OH) satisfying Eqs. p^T]) . fFM . 



region regardless to other parameters. 

First, let us consider t e . We expect that r E < A or t e > ^4 is the criterion for 
the scaling regime given by (|2-2ip and (|2-22p . where A is a constant. However, since 
r E is not a monotonic function of the area fraction eft and the restitution coefficient 
e, as shown in Fig. fT9t a). we conclude that neither r E < A or Tg > A is appropriate 
for the criterion. 

Similar to the case of r^., r* E in not a monotonic function of eft and e, as shown 
in Fig. rrW b). Therefore, we conclude that t* e is not an appropriate dimensionless 
time for the criterion. 
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Fig. 19. (a) 4> dependence on r E and (b) <j> dependence on t* e , for various e, from simulations 
with t* = 1.1 x 1CT 5 . 



Finally, let us consider r* and r* w , which are respectively related with r E and 
t* e as t* as 2r|/(l — e 2 ) and r* w « (1 — e 2 )r* E /2 in the collisional regime, but their 
dependency on (ft and e differs from those of t* e and r* s , as shown in Figs. l2"0T a) 
and [20T b) . Both r* and r*^ are monotonic functions of (j> and e. Since Eqs. (|2-21j) 
and (|2-22[) are satisfied in the high density region and r* and r* w are respectively 



Divergence of pressure and viscosity ... for hard and soft granular materials 21 



decreasing and increasing functions of the density cb, r* < A and > A are the 
possible conditions for the scaling given by Eqs. (|2-2ip and f|2-22j) . These conditions 
are also consistent with the dependencies of t* and t* u on e. Indeed, r* increases as 
the restitution constant increases, and r* w is a decreasing function of e. This means 
that the regions satisfying t* < A and r* w > A are narrower as the restitution 
constant increases, which is consistent with the numerical observation. Therefore, 
r * < A and r* w > A are the only two possible candidates to characterize the system 
with respect to their scaling behavior. It should be noted that r* w tends to zero in 
the hard disk limit r* — > 0. In this sense, to use r* w might involve a conceptual 
difficulty, even though r* w is finite in the jamming region. 
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Fig. 20. (a) <j) dependence on r* and (b) 

T* = 1.1 X 10" 5 . 



dependence on t* w for various e, from simulations with 



§4. Conclusion and Discussion 

In conclusion, we have investigated the dimensionless pressure P* and the di- 
mensionless viscosity n* of two-dimensional soft disk systems and have payed special 
attention to the rigid-disk limit of inelastically interacting systems, while near-rigid 
disks still have some elasticity ("softness"). 

For monodisperse systems, as the system approaches the elastic limit, e — > 1, 
both P* and rj* for cb < <f)„ = 0.71 approach the results of elastic rigid-disk systems, 
where the viscosity increases rapidly around cb = cb n due to ordering (crystallization) 
effects, while the pressure for cb > cb^ is still finite. 361 This result is consistent with 
Ref. l52|) . where Mitarai and Nakanishi suggested that the behavior of soft-disks in 
dilute collisional flow converges to that of rigid-disks in the rigid-disk limit. 

For polydisperse systems, both P* and rf behave as (cbj — cb)~ 2 and (cbj — <ft)~ 3 
near the jamming transition point, cbj > cj> v , as predicted in Refs. I39|) J4"0~|) . However, 
as the restitution coefficient e approaches unity, the scaling regime becomes narrower, 
and the exponents for the divergence of P* and rf' approach values close to —1 in 
the almost elastic case. 

From these results, we conclude that the predictions for the inelastic soft-disk 
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systems in Refs. I39]). l40p are applicable to the inelastic near-rigid disk systems below 
the jamming transition point, but the prediction cannot be used for almost elastic 
rigid-disk systems. It seems that r* w and r* are the only two possible candidates to 
characterize the criterion of this crossover. In other words, the energy dissipation rate 
and the shear rate set the two competing time-scales that define the dimensionless 
number r*. For r* <C 0.01 the near-rigid, dissipative scaling regime occurs, while for 
r* S> 0.01 the rigid, elastic scaling regime is realized. 

In three-dimensional sheared inelastic soft-sphere systems,®'^ even in monodis- 
perse cases, there is no indication of the strong ordering transition, and the scaling 
given in Eqs. (|2-2ip and (|2-22p seems to be valid. However, a direct comparison of 
near-rigid sphere with rigid sphere simulations in the spirit of the present study is 
unavailable to our knowledge. 

We restricted our interest to frictionless particles. When the particles have 
friction, the scaling relations for the divergence of the viscosity and the pressure 
may be different, as will be discussed elsewhere. Furthermore, the very soft or high 
shear rate regime also needs further attention in both 2D and 3D. 
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Appendix A 

The relation between Z and t* e 

In this appendix, we derive the relation between Z and r* E as 

z ~ t c r E l = T * cE , (A-i) 

which corresponds to the difference between counting contacts vs. counting of colli- 
sions in the simulations. (Note that counting contacts is not possible for rigid disks, 
since the probability to observe a t c = contact at any given snapshot in time is 
zero.) 

Since the ensemble average in Eq. (|3Tip is independent of i and j, without loss 
of generality, one can set i = 1 and j = 2, and obtains 

EE( K' " r ^)) = N(N - l)(0(a 12 - r 12 )>. (A-2) 

i j^i 
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Substituting this equation into Eq. (|3-lip . we obtain 

Z = (N -l)(0(a 12 -r 12 )). (A-3) 
On the other hand, t^ 1 is defined as the frequency of collisions per particle: 



1 



^ 1 T^oo T 



*E 1= E J im ^n c ,ij{T)\, (A-4) 



where n c ^j(t) is the number of the collisions between grains i and j until time T. 
Since limy_ i , 00 n c ,ij(T)/T is independent of j, like above, we obtain 

t E l = (N - 1) llim in«,i 2 (r) j (A-5) 

In order to derive Eq. (|A-lj) . the ensemble average in Eq. ()A-3p is replaced by 
the time average as 

Z = (N - 1) lim - [ T dt 0(a 12 - r 12 (t)), (A-6) 

where rij(t) is the distance between grains i and j at time t. Since 0(o~i 2 — ri 2 (t)) = 1 
for the duration t c after a collision begins, the integral in Eq. (|A-6|) is estimated as 
n c i2(T) t c , which yields 



Z = {N-l) lim l{n c , 12 (T)t c } 

T— »oo J 

(iV-l){ lim ^n C)12 (T) 



t c - (A-7) 



Finally, substituting Eq. (|A-5p into this equation, gives Eq. (|A- 1[) so that we can 
apply Eq. (|3 7 TTD . 
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